Time evolution of Matrix Product States 
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In this work we develop several new simulation algorithms for ID many-body quantum mechanical 
systems combining the Matrix Product State variational ansatz with Taylor, Pade and Arnoldi 
approximations to the evolution operator. By comparing with previous techniques based on Trotter 
decompositions we demonstrate that the Arnoldi method is the best one, reaching extremely good 
accuracy with moderate resources. Finally we apply this algorithm to studying how correlations 
are transferred from the atomic to the molecular cloud when crossing a Feschbach resonance with 
two-species hard-core bosons in a ID optical lattice. 



PACS numbers: 75.40.Mg, 02.70.-c, 03.75.Lm 

The Density Matrix Renormalization Group (DMRG) 
method is a successful technique for simulating large low- 
dimensional quantum mechanical systems . Developed 
for computing ground states of ID Hamiltonians, it is 
equivalent to a variational ansatz known as Matrix Prod- 
uct States (MPS) % This relation has been recently 
exploited to develop a much wider family of algorithms 
for simulating quantum systems .including time evolu- 
tion 13, Q , finite temperature 0, Q and excitation spec- 
tra p. Some of these algorithms have been translated 
back to the DMRG language using optimizations de- 
veloped in that field and introducing other techniques 
such as Runge-Kutta or Lanczos approximations of the 
evolution operator 0, • 

The MPS are a hierarchy of variational spaces, Sd, 
[See Eq. ^] sorted by the size of its matrices, D. MPS 
can efhciently represent many-body states of ID systems, 
even when the Hilbert space is so big that the coefficients 
of a pure state on an arbitrary basis cannot be stored in 
any computer. While the accuracy of this representa- 
tion has been proven for ground states IQJ , evolution of 
an arbitrary state changes the entanglement among its 
parties, and a MPS description with moderate resources 
(small D) might stop to be feasible. 

We will take a pragmatic approach. First of all, MPS 
and DMRG algorithms can compute truncation errors 
so that the accuracy of simulations remains controlled. 
Second, we are interested in simulating physically small 
problems, such as the dynamics of atoms and molecules 
in optical lattices. For such problems small D are suffi- 
cient to get a qualitative and even quantitatively good 
description of the observables in our systems. As we 
will see below, the biggest problem is the accumulation 
of truncation errors and not the potential accuracy of a 
given space So for representing our states. 

The outline of the paper is as follows. We first present 
many recently developed simulation algorithms under a 
common formalism based on the optimal truncation op- 
erator. We then introduce two new algorithms: one of 
them is based on Pade expansions of the evolution op- 
erator while the other one uses linear combinations of 
MPS to increase the accuracy (similar to an "Arnoldi" 
method). A comparison of all algorithms using spin— ^ 



models shows that all methods are strongly limited by 
truncation and rounding errors and that the Arnoldi 
method performs best for a given accuracy. 

In the last part of the paper we apply these algo- 
rithms to study a model of hard-core bosonic atoms go- 
ing through a Feschbach resonance. Current experiments 
1 11] with such systems have focused on the number and 
stability of the formed molecules. In this work we focus 
on the ID many-body states and show that coherence is 
transferred from the atomic component to the molecular 
one, so that this procedure can be used to probe higher 
order correlations in the atomic cloud. 

Summary of algorithms. -"Die space of MPS of size D is 
the set of states of the form 



Sn ■■= 
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(1) 



Here ik ~ 1 . . .d labels the physical state of the A:-th 
lattice site, the Hilbert space has a total of states 
and we sum over repeated indices. 

Since So is not a vector space, we cannot solve a 
Schrodinger equation directly on it. Our goal is rather to 
approximate the evolution at short times by a formula of 
the kind + Ai)) Vd [Un{M)\ip{t))] . The operator 
Vd optimally projects an arbitrary vector onto the man- 
ifold Sd, and C/„(At) = exp{-iHAt) + 0{Ar') is itself 
an approximation to the evolution operator. 

This formulation applies qualitatively to all recently 
developed MPS and DMRG algorithms. For instance, 
Refs. 3, 4, "Tl use a Trotter decomposition of a nearest 
neighbor Hamiltonian H = J2m Hm.m-i-i 



U2{At) = e 



2At 



xe 



-iEfc -H2fc,2fc+lAt/2 



(2) 



while the Runge-Kutta algorithm in Ref. ^ is equiv- 
alent to a fourth order expansion of the exponential 
U4{At) — J2ti=o (~*^^^)"- practice there are dif- 
ferences between all cited implementations, such as trun- 
cating in between applications of the exponentials or re- 
placing the previous Taylor expansion by a product of 
truncated first order monomials, Y[k ^-d(1 + ctkH). 
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FIG. 1: (a) Error, e, vs. time step, At, for simulations 
with 8 spins, D = 16 and T = 10. We compare Trotter 
expansions (circles), Runge-Kutta of 4th order (red). Fade 
expansions of 2nd and 4th order (blue) and Arnoldi with 4 
and 8 vectors (black). In dash-dot line we plot the estimate 
£ = (r/At)^10"^'^. (b) Same as plot but for 16 spins and 
comparing only Arnoldi-4, Trotter and Arnoldi-8 (top to bot- 
tom) for D = 64,80 and 128 (dash, dash-dot, solid). While 
(a-b) are done for Q := arctan(i3/J) = 7r/4, A = 0, we have 
scanned the whole Hamiltonian space, 6 £ [0, 27r], A G [—3, 3]. 
In (c-d) we show sections along d and along A, using 16 spins, 
T = 5, -D = 32 and the methods from (b). 



In this work we present two new methods. The first 
method uses a Pade approximation to the exponential 



C/n(Ai) 



(3) 



computed with same order polynomials in the numerator 
and denominator [T^. 

The second method combines ideas from Refs. [^|^. 
First of all, it is possible to prove that a linear combi- 
nation of MPS, such as Cfc|0^''), resides in a space 
of bigger matrix product states, Sn^d- In the language 
of DMRG, vectors each of size D provide us with an 
effective basis of size N^D"^ . This optimistic estimate is 
only possible when the vectors |<^^^) are linearly inde- 
pendent. Our choice will be 



k+i)c^rD\H\c^k)-Y. 

3<k 



{4'j\H\M 



10.) 



(4) 



with initial condition \(j>o) |'0(O))- Finally, defining the 
matrices Nih := {4>j\4>i) and Hik '■— {4>j\H\(j)i) we com- 
pute an Arnoldi estimate of the exponential 



(5) 



This algorithm involves several types of errors. The er- 
ror due to using only Ny basis vectors is proportional the 
norm of the vector as in Lanczos algorithms 
Truncation errors arising from Vd can be computed and 
controlled during the numerical simulation. It is impor- 
tant to remark that while the errors in Eq. |0J may be 
compensated by adding more vectors, the biggest errors 
arise from the final truncation in Eq. 

Implementation.- All algorithms build on the same 
nonlinear operator, Vd, which optimally approximates 
a state using the elements of Sd- To truncate a linear 
combination of vectors we use the definition IJl 



argmm 



(6) 



If we rather want to approximate the action of an opera- 
tor that can be decomposed as J7 = X^^Y, we will apply 
a generalization of the correction vector method |^ 

Vd {X-'Y\<f>)) := argmin||X|7^) - y|0)||. (7) 

This last formula has been used for the Pade J^Jl equation 
and for the Runge-Kutta methods, where X and Y are 
polynomials of the Hamiltonian. 

One may quickly devise a procedure for computing the 
minimum of Eq. Q based on the definition of distance: 



\\X\^)-Y 



{tl;\X^X\tP)-2^{tP\X''Y[ 



\Y\ 



(8) 

All scalar products in Eq. © are quadratic forms with 
respect to each of the matrices in the states |(/)) and 1-0). 
The distance is minimized by optimizing these quadratic 
forms site by site, or two sites at a time, sweeping over 
all lattice sites until convergence to a small value which 
will be the truncation error 

The performance our algorithms is upper bounded by 
0{NyNHD^ / L), where the Nh is related to the number 
of operators in {X, Y, X'^X, Y'^Y, . . .}, L is the size of the 
problems and Ny is the number of vectors in Eq. © . For 
open boundary conditions the cost reduces by 0{D^). 

Comparison.- We tested all algorithms by simulating 
the evolution under the family of spin-i Hamiltonians 
with nearest-neighbor interactions 
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Aslsl+,) + Bsl] , (9) 



of an initial state, |V'(0)) oc ([0) -\- , where |0) and 

1) are the eigenstates of s^. By restricting ourselves to 
"small" problems {L < 16), we could compare all algo- 
rithms with accurate solutions based on exact diagonal- 
izations and the Lanczos algorithm p^ . 

In Fig. ^ we plot the errors of the different meth- 
ods, e := Wi^oiT) - C/(T)V'(0)1P, for 8 spins, D = 16, 
J = 2B and A = 0. Since So contains all possible states, 
Vd = 1 and we expect no truncation errors. There- 
fore, for medium to long time steps the Trotter, Taylor 
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FIG. 2: (a) Error made by approximating exp(— i_fff)|i/)(0)) 
with a vector in 5*64 (solid) and error made when simulating 
this Schrodinger equation using the Arnoldi method with 8 
vectors, D = QA and /St — 0.16 (dashed), (b) Similar as 
before, errors in the Arnoldi method for varying number of 
vectors, T = 10 and At = 0.04, 0.16, 0.32 (bottom to top). 



and Pade approximations show the expected behavior, 
with an error growing like 0{At^) or ©(At^). The er- 
rors of the Arnoldi decompositions qualitatively follow 
0{At^^). However, all methods break these ideal laws at 
some point, acquiring an error of order 0{At^^), which 
is exponential in the number of steps T/ At. This growth 
signals the finite accuracy of the optimization algorithms, 
due to the limited precision of the computer. Roughly, 
since current computers cannot compute the norms of 
vectors, with a relative error better than 10~^^, a 

worst case estimate is e = {T/AtflQ-^^ [See Fig.[Bi]. 

We have compared the best methods, that is the Trot- 
ter decomposition and the Arnoldi methods with 4 and 8 
vectors, using a bigger lattice with 16 spins, and smaller 
matrices, D = 64, 80 and 128. As Fig. shows, the 
XY model establishes entanglement along the lattice, the 
smaller matrices cannot account for this and the errors 
grow exponentially. Figure|3i. shows how the errors in the 
Arnoldi method are correlated to the errors made when 
approximating the exact solution with a MPS of fixed 
size. One could think that by increasing the number of 
vectors we should be able to also increase the integration 
time-step and decrease the truncations even for fixed D, 
but as Fig. |3) shows, this not entirely true. 

Summing up, one should use the method that allows 
for the longest time steps and the least number of trun- 
cations (or applications of Vd) and mathematical oper- 
ations. All methods have an optimal time step which is 
a compromise between the errors in J7„ and the round- 
ing and truncation errors made on each step. Regarding 
performance, the Arnoldi decomposition is competitive 
with the Trotter decomposition as it reaches the same 
accuracy for longer time steps. Furthermore, there is a 
huge potential for parallelizability, roughly 0{NyNH / L), 
which all other presented algorithms lack. 

The fact that previous results are model-independent is 
confirmed by a systematic scanning of all possible Hamil- 
tonians in Eq. Q . A selection is shown in Fig. nj;-d. The 
Arnoldi method is shown to be the most accurate one, 
even for gapless problems. When the Arnoldi method 
fails it is due to truncation errors — the evolved state 
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FIG. 3: Dynamics of a site with two atoms when the energy of 
the molecules is ramped linearly: Em = f/|| -I- 4f2(l — 2t/T). 
We plot (a) the instantaneous energy levels (solid), and (b) 
the fraction of atoms converted into molecules. 



cannot be accurately represented by MPS — , which af- 
fect equally aU MPS/DMRG algorithms [Fig. [Hi]. 

Molecules in the lattice - We have used the Trotter and 
Arnoldi methods to simulate the conversion of bosonic 
atoms into molecules, when confined in an optical lattices 
and moving through a Feschbach resonance |l3|. The 
goal is to study how correlation properties are transferred 
from the atoms into the molecules and how this dynamics 
is affected by atom motion and conversion efficiency. 

The effective model combines the soft-core Bose- 
Hubbard model used to describe the Tonks gas exper- 
iments 14] , with a coupling to a molecular state |15| 



H = -J a\^aj„ + ^ ^^^a\^a\^a^„a^„ (10) 



^ [{E„, + C/„n|"))r4™^ + r![6la»T«a + H. c.]| 



Here, a^i, Oij^ and bi are bosonic operators for atoms in 

two internal states and the molecule; rS"'^ and nf^^ are 
the total number of atoms and of molecules on each site, 
and we have the usual two-level coupling with Rabi fre- 
quency Q. and detuning A :— E„i — Um- For simplicity, we 
will assume that atoms and molecules interact strongly 
among themselves (?7f,T' oo), so that we can 

treat them as hard-core, of ^,bf,airjbi — 0. Also since 
molecules are heavier, we have neglected their tunneling 
amplitude, although that could be easily included. 

As the energy of the molecular state is shifted from 
E,n » Um down to Em <^ Um, the ground state of 
Eq. H10|l changes from a pair of coupled of Tonks gases, 
to a purely molecular insulator. We want to study the 
dynamics of this crossover as E„i is ramped slowly from 
one phase to the other. 

The simplest situation corresponds to no hopping: iso- 
lated atoms experience no dynamics, while sites with 
two atoms may produce a molecule. The molecular and 
atomic correlations at the end of the process are directly 
related to two-body correlations in the initial state |16j . 



{m\mk)t=T 
{alak)t=T 



(11) 



{alak)t=o - {nk^nki)t=o- (12) 
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FIG. 4: (a) Fraction of atoms converted into molecules vs. 
adimensionalized hopping amplitude, for U^^i = and U^i — 
1, from top to bottom. Circles show the outcome of the nu- 
merical experiment, while solid lines contain the ideal fraction 
lim . The case J — uses an initial condition J = O.lfi and 
then switches off tunneling before ramping, (b) Correlations 
of the molecular state, {m\mj)t=T (circles) after the ramp, 
and those of the initial atomic state, (al|fii|ffljT'*ii)*=o (solid 
line). The plot for = 5 J has been shifted up by 0.2. 

Therefore, this process can thus be used as a tool to probe 
quantum correlations between atoms. Studying the two- 
level system {a[.f aj^jO), 6^|0)}, we conclude that for this 
process to work with a 90% efficiency, the ramping time 
should be larger than T ~ 1.5/n [See Fig.EIb)]. 

We have simulated numerically the ramping of small 
lattices, L = 10 to 32 sites, with an initial number of 
atoms N-f^i — L/2, 3L/4. The value of the molecular cou- 
pling has been fixed to = 1 and the interaction has been 
ramped according to _E,„ — U-\i + 457(1 — 2t/T) using the 
ideal ramp time T — 1.5/57. We have used two particular 
values of the inter-species interaction, U^^ i/ft = 0, 2, and 
scanned different values of the hopping J/fl g [0,0.4]. 
The initial condition was always the ground state of the 
model with these values of Un/ J and no coupling. These 



states contain the correlations that we want to measure. 

The main conclusions is that indeed the correlations 
of the molecules are almost those of the initial state of 
the atoms Hll|l . even for J = 0.457 when the process 
has not been adiabatic. An intuitive explanation is that 
hopping is strongly suppressed as we approach the reso- 
nance, due to the mixing between atomic states, which 
can hop, and molecular states, which are slower. We can 
say that the molecules thus pin the atoms and measure 
them. This explanation is supported by a perturbation 
analysis at J <C 57, where one finds that a small molec- 
ular contamination slows the atoms on the lattice. This 
analysis breaks down, however, for J ~ 57, the regime in 
which the numerical simulations are required. 

Conclusions.- We have performed a rather exhaustive 
comparison of different methods for simulating the evolu- 
tion of big, one-dimensional quantum systems 
with two other methods developed in this work. All 
methods are substantially affected by truncation and 
rounding errors, and to fight the latter we must choose 
large integration time-steps. However, the only algorithm 
which succeeds for large time-steps is an Arnoldi method 
developed in this work. Finally, this algorithm can be 
applied to problems with long range interactions. 

Using this algorithm, we have simulated the dynam- 
ics of cold atoms in a ID optical lattice when crossing a 
Feschbach resonance. The main conclusion is that with 
rather fast ramp times it is possible to map the correla- 
tions of the atomic cloud (two Tonks gases in this case) 
and use this as a measuring tool in current experiments. 
This result connects with similar theoretical predictions 
for fermions in Ref. O. Simple generalizations of this 
work will allow us in the future to analyze losses and 
creation of entanglement with the help of the molecular 
component. 
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